OHI-Northeast | OHI Science | Citation policy
OHI-Northeast | OHI Science | Citation policy
Downloaded: December 14, 2018 (emailed to us by Jeffrey Vieser at NMFS)
Description: Records of Bmsy and Fmsy estimates from stock assessments conducted in the greater Northeast region
Time range: 2004 - 2018
Format: Tabular
Load raw data
#stock assessment data shared by NMFS
#northeast region
ne_data <- readxl::read_excel(file.path(dir_anx, '_raw_data/NOAA_NMFS/stock_assessments/NE_Assessment_Records.xlsx'))
#midatlantic region
midatl_data <- readxl::read_excel(file.path(dir_anx, '_raw_data/NOAA_NMFS/stock_assessments/MA_Assessment_Records.xlsx'))Clean data
#select just the columns we are interested in
data <- ne_data %>%
bind_rows(midatl_data) %>%
select(year = Year,
stock = Stock,
f_fmsy = `F/Fmsy`,
b_bmsy = `B/Bmsy`) %>%
gather(key = "indicator", value = "value", -year, -stock) %>%
filter(value > 0) #Smooth skate has F/Fmsy value of -9999 so removing that
ggplot(data, aes(x = year, y = value, color = indicator)) +
geom_hline(yintercept = 1, color = "darkgray") +
geom_line() +
facet_wrap(~stock, labeller = labeller(stock = label_wrap_gen(width = 25)), scales = "free_y") +
theme_bw() +
theme(strip.text = element_text(size=8))
Save data
write.csv(data, file = "data/nmfs_stock_assessment_data.csv")Span plot
stock_ids <- data$stock %>% unique()
spanplot_df <- data %>%
group_by(stock) %>%
mutate(yr_min = min(year), yr_max = max(year))
span_plot <- ggplot(spanplot_df, aes(x = stock)) +
#ggtheme_basic +
theme(panel.border = element_blank(),
panel.grid.minor = element_blank(),
panel.grid.major = element_line(colour = 'grey90'),
panel.background = element_blank()) +
geom_linerange(aes(ymin = yr_min, ymax = yr_max,
color = indicator),
position = position_dodge(width = .5), alpha = .8) +
scale_color_manual(values = c('darkred', 'darkblue'),
labels=c("B/Bmsy", "F/Fmsy")) +
coord_flip() +
labs(x = 'Stock',
y = 'Year',
color = 'Indicator')
span_plot
What species have multiple stocks? These ones will give us a tough time when matching up catch to the stocks.
dups <- data %>%
separate(stock, into = c("species", "location"), sep = " - ") %>%
select(-indicator, -value) %>%
distinct() %>%
group_by(year, species) %>%
mutate(count = n()) %>%
ungroup() %>%
filter(count > 1) %>%
mutate(stock = paste0(species, " - ", location)) %>%
select(stock, species, location) %>%
distinct()
dups## # A tibble: 21 x 3
## stock species location
## <chr> <chr> <chr>
## 1 Atlantic cod - Georges Bank Atlantic cod Georges Bank
## 2 Atlantic cod - Gulf of Maine Atlantic cod Gulf of Maine
## 3 Haddock - Georges Bank Haddock Georges Bank
## 4 Haddock - Gulf of Maine Haddock Gulf of Maine
## 5 Windowpane - Gulf of Maine / Geo… Windowpane Gulf of Maine / George…
## 6 Windowpane - Southern New Englan… Windowpane Southern New England /…
## 7 Winter flounder - Gulf of Maine Winter flound… Gulf of Maine
## 8 Winter flounder - Southern New E… Winter flound… Southern New England /…
## 9 Yellowtail flounder - Cape Cod /… Yellowtail fl… Cape Cod / Gulf of Mai…
## 10 Yellowtail flounder - Georges Ba… Yellowtail fl… Georges Bank
## # … with 11 more rowsLooks like there are just 8.
Use a linear model to fill years between assessments.
We gapfill at B/bmsy not stock score. If we have 2010 and 2012 data that shows b/bmsy of 0.7 and 1.3 respectively. If we gapfill at B/bmsy, the 2011 value becomes 1. This would get a stock score of it. But if we instead gapfill between 2010 anda 2012 stock scores (0.7 and something like 0.9) we would get 0.8…
data_gf <- data %>%
group_by(stock, indicator) %>%
mutate(timeseries = n()) %>% #get number of years for each stock/indicator combo. We need at least 2 years
filter(timeseries > 1) %>%
complete(year = 2005:2018) %>%
ungroup() %>%
mutate(value2 = ifelse(year == 2018 & is.na(value), zoo::na.locf(value), value), #if 2018 is NA, use 2017 (or last reported value)
value3 = ifelse(year == 2005 & is.na(value2), zoo::na.locf(value2, fromLast = TRUE), value2)) %>% #if 2005 is NA, roll back the closest value from following years
group_by(stock, indicator) %>%
mutate(gf_value = zoo::na.approx(value3),
gf = ifelse(is.na(value), 1, 0)) %>% #if a value is gapfilled, assign it a 1, otherwise 0
select(-value, -value2, -value3, -timeseries) %>%
rename(value = gf_value)Convert these value to stock scores.
#grab just B/Bmsy data
nmfs_b_bmsy <- data_gf %>%
filter(indicator == "b_bmsy") %>%
select(stock, year, value)
#grab just F/Fmsy data
nmfs_f_fmsy <- data_gf %>%
filter(indicator == "f_fmsy") %>%
select(stock, year, value)
###
b_bmsy_underexploit_penalty <- 0.25
b_bmsy_underexploit_thresh <- 3.00
f_fmsy_underfishing_penalty <- 0.25
f_fmsy_overfishing_thresh <- 2.00
### Apply rolling mean to F/Fmsy
## Why do we do this? because B is a less sensitive metric (relies of biological processes) and F can fluctuate pretty easily because it is really just a mgmt decision.
nmfs_f_fmsy <- nmfs_f_fmsy %>%
arrange(stock, year) %>%
group_by(stock) %>%
mutate(f_fmsy_rollmean = zoo::rollmean(value, k = 3, align = 'right', fill = NA)) %>%
ungroup() %>%
select(-value) %>%
rename(value = f_fmsy_rollmean)
stock_status_layers <- nmfs_b_bmsy %>%
full_join(nmfs_f_fmsy) %>%
spread(indicator, value)
########################################################.
##### run each fishery through the Kobe plot calcs #####
########################################################.
### * ram_b_bmsy, ram_f_fmsy
### Function for converting B/Bmsy values into a 0 - 1 score
rescale_bprime_crit <- function(fish_stat_df,
bmax, bmax_val) {
###using NOAA's limits here
overfished_th <- 0.8
###
underfished_th <- 1.2
bmax_adj <- (bmax - underfished_th) / (1 - bmax_val) + underfished_th
### this is used to create a 'virtual' B/Bmsy max where score drops
### to zero. If bmax_val == 0, this is bmax; if bmax_val > 0, bmax_adj
### extends beyond bmax, to create a gradient where bmax_val occurs at bmax.
fish_stat_df <- fish_stat_df %>%
# group_by(stock) %>% ### grouping by stock will set b_max by max per stock, instead of max overall
mutate(b_max = max(b_bmsy, na.rm = TRUE)) %>%
ungroup() %>%
mutate(bPrime = NA,
bPrime = ifelse(b_bmsy < overfished_th, ### overfished stock
b_bmsy / overfished_th,
bPrime),
bPrime = ifelse(b_bmsy >= overfished_th & b_bmsy < underfished_th,
1, ### appropriately fished stock
bPrime),
bPrime = ifelse(b_bmsy >= underfished_th,
(bmax_adj - b_bmsy) / (bmax_adj - underfished_th), ### underfished stock
bPrime),
bPrime = ifelse(bPrime < 0, 0, bPrime))
return(fish_stat_df)
}
### Function to create vertical gradient based on distance from
### ideal F/Fmsy value to actual F/Fmsy value
f_gradient <- function(f, over_f, under_f, fmax, fmin_val) {
x <- ifelse(f < over_f & f > under_f, 1, NA)
x <- ifelse(f <= under_f, (f * (1 - fmin_val) / under_f + fmin_val), x)
x <- ifelse(f >= over_f, (fmax - f) / (fmax - over_f), x)
x <- ifelse(f > fmax, NA, x)
return(x)
}
### Function to convert F/Fmsy values into 0 - 1 score
rescale_fprime_crit <- function(fish_stat_df,
fmax, fmin_val) {
### params - taken from BC but changed Bcrit to 0 instead of 0.4:
Bcrit <- 0.5; overfished_th <- 0.8
### underfishing_th is set to the idea "1/3 for the birds":
underfishing_th <- 0.66; overfishing_th <- 1.2
bcritslope = 1 / (overfished_th - Bcrit)
### connecting from (Bcrit, 0) to (overfished_th, 1)
fish_stat_df <- fish_stat_df %>%
mutate(fPrime = ifelse(b_bmsy < overfished_th & f_fmsy < fmax,
f_gradient(f_fmsy + (overfished_th - b_bmsy) * bcritslope,
over_f = overfishing_th,
under_f = underfishing_th,
fmax = fmax,
fmin_val = fmin_val),
NA),
fPrime = ifelse(b_bmsy >= overfished_th & f_fmsy < fmax,
f_gradient(f_fmsy,
over_f = overfishing_th,
under_f = underfishing_th,
fmax = fmax,
fmin_val = fmin_val),
fPrime),
fPrime = ifelse(is.na(fPrime), 0, fPrime), ### fill zeros everywhere unscored
fPrime = ifelse(is.na(f_fmsy), NA, fPrime) ### but if no f_fmsy, reset to NA
)
return(fish_stat_df)
}
stock_status_df <- stock_status_layers %>%
rescale_bprime_crit(bmax = b_bmsy_underexploit_thresh,
bmax_val = b_bmsy_underexploit_penalty) %>%
rescale_fprime_crit(fmax = f_fmsy_overfishing_thresh,
fmin_val = f_fmsy_underfishing_penalty) %>%
mutate(x_prod = ifelse(!is.na(fPrime), (fPrime * bPrime), bPrime),
basis = case_when(
!is.na(fPrime) & !is.na(bPrime) ~ 'F_Fmsy, B_Bmsy',
is.na(fPrime) & !is.na(bPrime) ~ 'B_Bmsy only',
is.na(bPrime) & !is.na(fPrime) ~ 'F_Fmsy only'
)) %>%
dplyr::select(year, stock,
score = x_prod,
basis,
bPrime, fPrime,
b_bmsy, f_fmsy) Take the average scores for the 5 species with sub-stocks. Just doing a group_by with species and averaging scores will do this for all species (but only change the values for these sub-species).
Remove stocks with just F_Fmsy - these stocks have no stock scores since we don’t have a way to get a score from just F/Fmsy. Also roll all values forward to 2017.
stock_scores <- stock_status_df %>%
filter(basis != "F_Fmsy only") %>%
group_by(stock, year) %>%
mutate(score = mean(score, na.rm = T)) %>%
select(year, stock, score, bPrime, fPrime, b_bmsy, f_fmsy, basis)
write.csv(stock_scores, file = "data/stock_scores_nmfs.csv") ggplot(stock_scores, aes(x = year, y = score)) +
geom_line() +
theme_bw() +
facet_wrap(~stock)
kobe_plot <- function(sp){
ss_df <- stock_scores %>%
filter(stock == sp) %>%
arrange(year) %>%
mutate(last_bbmsy = last(b_bmsy),
last_ffmsy = last(f_fmsy),
last_datayear = last(year)) %>%
ungroup()
generate_kobe_df <- function(f_fmsy_max = 2.5,
b_bmsy_max = 3.0,
reso = 0.01,
bmax_val = 0,
fmin_val = 0,
weighting_b = 1) {
kobe_raw <- data.frame(stock = 1,
f_fmsy = rep(seq(0, f_fmsy_max, reso), each = round(b_bmsy_max/reso) + 1),
b_bmsy = rep(seq(0, b_bmsy_max, reso), times = round(f_fmsy_max/reso) + 1))
kobe <- kobe_raw %>%
rescale_bprime_crit(bmax = b_bmsy_underexploit_thresh,
bmax_val = bmax_val) %>%
rescale_fprime_crit(fmax = f_fmsy_overfishing_thresh,
fmin_val = fmin_val) %>%
mutate(x_geom = (fPrime * bPrime),
x_arith = (fPrime + bPrime) / 2)
return(kobe)
}
bbmsy_lim <- max(round(max(ss_df$b_bmsy, na.rm = TRUE) + .1, 1), 3)
ffmsy_lim <- max(round(max(ss_df$f_fmsy, na.rm = TRUE) + .1, 1), 2.5)
kobe_df <- generate_kobe_df(f_fmsy_max = ffmsy_lim,
b_bmsy_max = bbmsy_lim,
bmax_val = .25,
fmin_val = .25)
kobe_stock_plot <- ggplot(data = kobe_df, aes(x = b_bmsy, y = f_fmsy)) +
theme_bw() +
geom_raster(alpha = .8, aes(fill = x_geom)) +
scale_fill_distiller(palette = 'RdYlGn', direction = 1) +
labs(title = as.character(sp),
x = 'B/Bmsy',
y = 'F/Fmsy',
fill = "Stock score") +
geom_path(data = ss_df,
show.legend = FALSE,
aes(x = b_bmsy, y = f_fmsy, group = sp),
color = 'grey30') +
geom_point(data = ss_df,
show.legend = FALSE,
aes(x = last_bbmsy, y = last_ffmsy)) +
geom_text(data = ss_df %>%
mutate(year = ifelse(year/5 == round(year/5) | year == last_datayear, year, NA)),
aes(x = b_bmsy, y = f_fmsy, label = year),
hjust = 0, nudge_x = .05, size = 2)
return(kobe_stock_plot)
}Apply function to all species in the stock_scores data frame.
sp <- stock_scores %>%
filter(basis == "F_Fmsy, B_Bmsy")
plots <- lapply(unique(sp$stock), kobe_plot)
cowplot::plot_grid(plotlist = plots[1:4])
cowplot::plot_grid(plotlist = plots[5:8])
cowplot::plot_grid(plotlist = plots[9:12])
cowplot::plot_grid(plotlist = plots[13:16])
cowplot::plot_grid(plotlist = plots[17:20])
cowplot::plot_grid(plotlist = plots[21:24])
cowplot::plot_grid(plotlist = plots[25:28])
cowplot::plot_grid(plotlist = plots[29:32])
cowplot::plot_grid(plotlist = plots[33:35])
write.csv(stock_scores, file = "data/nmfs_stock_scores.csv")